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Abstract 

An  orbit  that  lies  on  a  Kolmogorov,  Arnold,  and  Moser  (KAM)  Torus  will 
remain  on  that  torus  until  and  unless  it  experiences  a  force  that  causes  it  to  leave 
the  torus.  Earth  satellites  that  are  subject  only  to  the  Earth’s  gravity  held  may 
lie  on  such  KAM  tori.  Analyzing  on  orbit  satellite  position  data  should  allow  for 
the  identification  of  the  fundamental  frequencies  needed  to  define  the  KAM  tori  for 
modeling  Earth  satellite  orbits. 

KAM  Tori  are  created  for  the  Gravity  Recovery  and  Climate  Experience  (GRACE) 
and  Jason-1  satellites  to  model  their  orbital  motion.  Precise  position  data  for  the 
satellites  is  analyzed  using  a  modified  Laskar  frequency  algorithm  to  determine  the 
fundamental  frequencies  of  the  orbits.  The  fundamental  frequencies  along  with  a  set 
of  Fourier  coefficients  completely  describe  the  tori.  These  tori  are  then  compared  to 
the  precise  orbital  position  data  for  the  satellites  to  determine  how  well  they  model 
the  orbits. 

The  KAM  torus  created  for  the  Jason-1  satellite  is  able  to  represent  the  position 
of  the  satellites  to  within  1  km.  Further  refinement  of  the  torus  should  be  possible, 
resulting  in  a  more  accurate  model  of  the  orbit.  The  GRACE  torus  was  less  successful 
at  determining  the  satellite  positions.  Atmospheric  drag  cannot  be  ignored  at  the 
altitude  where  GRACE  flies.  It  may  still  be  possible  to  model  GRACE  with  a  KAM 
torus  by  applying  perturbation  theory  to  the  torus;  however,  further  research  is  needed 
to  confirm  this. 
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Application  of  KAM  Theorem  to  Earth  Orbiting 

Satellites 


I.  Introduction 

Since  the  beginning  of  the  space  age,  thousands  of  objects  have  been  launched 
into  orbit.  However,  unlike  terrestrial  and  airborne  systems  that  are  moved  to  places 
where  they  are  not  going  to  interfere  with  ongoing  operations,  space  systems  continue 
to  float  in  space  until  their  orbit  decays  enough  to  reenter  the  atmosphere.  Thus, 
there  are  still  many  non-operational  systems  flying  in  space,  along  with  our  operational 
systems.  The  US  Space  Surveillance  Network  (SSN)  tracks  these  objects  and  attempts 
to  accurately  predict  their  future  motion  to  ensure  that  high  priority  missions  are  not 
impacted. 

1 . 1  Motivation 

The  standard  method  for  orbit  modeling  starts  with  Kepler’s  solution  to  the 
Two  Body  Problem  (2BP)  and  uses  perturbation  theory  to  account  for  some  of  the 
major  error  sources;  these  errors  are  due  to  effects  like  variations  in  the  Earth’s  grav¬ 
itational  field.  A  common  method  begins  with  an  initial  condition,  and  the  orbits 
are  integrated  forward  in  time.  This  is  a  computationally  expensive  process,  which 
restricts  the  number  of  perturbations  that  can  be  included;  for  example,  the  Earth’s 
oblateness  effect  is  typically  included,  but  other  terms  of  the  geopotential  may  be  left 
out,  especially  when  integrating  large  numbers  of  orbits.  Other  methods  are  analyti¬ 
cal  and  use  series  expansions  of  the  equations  of  motion,  which  include  perturbations 
like  the  geopotential.  These  series  expansions  become  very  complex  and  must  be  trun¬ 
cated  to  make  them  easier  to  handle  operationally.  Orbital  predictions  created  with 
these  methods  are  only  valid  for  a  short  time  before  they  must  be  recalculated;  the 
perturbations  that  are  not  accounted  for,  like  the  higher  order  gravitational  variations, 
build  up  in  the  error  over  time. 
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1.2  Approach 

A  more  accurate  method  for  orbit  determination  could  reduce  the  work  required 
to  track  all  of  the  objects  orbiting  Earth.  This  work  attempts  to  show  that  the  KAM 
theorem  may  be  one  such  method  that  could  be  used  to  model  satellite  orbits. 

The  Kolmogrov-Arnold-Moser  (KAM)  theorem  states  that  a  nearly  integrable 
Hamiltonian  system  subject  to  a  small  perturbation  will  lie  on  an  invariant  torus  [9]; 
these  tori  have  become  known  as  KAM  tori.  Earth  orbiting  satellites  may  behave 
according  to  the  KAM  theorem  since  the  2BP  is  an  integrable  Hamiltonian  system. 
Assuming  the  only  perturbation  source  is  the  Earth’s  gravitational  potential,  which 
causes  small,  smooth  perturbations  to  the  2BP  [19],  the  criteria  for  KAM  theorem  is 
met. 

The  key  to  determining  whether  a  system  lies  on  a  KAM  torus  lies  in  finding 
the  fundamental  frequencies  of  the  torus  [9];  the  number  of  frequencies  corresponds  to 
the  number  of  coordinates  in  the  Hamiltonian.  In  the  case  of  Earth  orbiting  satellites 
there  are  three  fundamental  frequencies:  the  anomalistic  frequency  (cni),  the  Earth 
rotation  rate  plus  the  nodal  regression  rate  (u;2)  and  the  apsidal  regression  rate  (u;3). 
These  frequencies  can  be  identified  from  a  Fourier  transform  of  the  orbital  position 
data  of  a  satellite.  If  the  satellite’s  orbit  lies  on  a  KAM  torus,  the  satellite  position 
can  be  determined  accurately  for  any  point  in  the  future,  assuming  no  other  external 
forces  are  encountered. 

1.3  Problem  Statement 

This  work  looks  at  two  satellite  systems,  the  Gravity  Recovery  and  Climate 
Experiment  (GRACE)  and  Jason-1  satellites.  Orbital  data  from  each  system  is  run 
through  a  modified  Laskar  Fourier  transform  algorithm  to  determine  the  fundamental 
frequencies  of  the  orbit  [10].  These  frequencies  and  their  corresponding  Fourier  series 
coefficients  define  the  new  coordinates  for  the  Hamiltonian  that  describes  the  torus; 
the  conjugate  momenta  can  be  found  from  the  coordinates  and  their  derivatives  [20]. 
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1.4  Results 

Analysis  shows  that  the  Jason-1  satellite  orbit  can  be  modeled  by  a  KAM  torus 
to  within  1  km  of  the  real  orbital  data  over  a  30  day  period.  Further  refinement  may 
improve  this  accuracy.  The  torus  for  the  GRACE  satellite  was  not  successful,  due 
to  the  appearance  of  atmospheric  drag  in  the  data.  It  may  still  be  possible  to  define 
the  GRACE  orbit  with  a  KAM  torus,  if  the  torus  can  be  modified  to  account  for  the 
atmospheric  drag. 
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II.  Background 

This  chapter  begins  with  discussions  of  the  GRACE  and  Jason-1  satellites.  These 
satellites  were  chosen  because  the  position  data  for  each  is  accurately  known.  Accurate 
position  data  is  needed  to  ensure  the  frequencies  can  be  identified  with  as  much 
accuracy  as  possible. 

Current  methods  for  orbit  determination  are  also  explained  along  with  their 
drawbacks.  The  KAM  theorem  will  be  discussed  in  more  detail,  and  past  applications 
are  presented.  Finally,  explanations  of  the  orbital  dynamics  and  the  geopotential  are 
given. 

2.1  GRACE 

The  mission  of  GRACE  is  to  measure  the  temporal  and  spatial  variability  of 
Earth’s  gravitational  field.  There  are  some  advantages  that  make  GRACE  supe¬ 
rior  to  previous  methods  for  creating  a  gravity  model:  the  use  of  identical  satellites 
flying  in  formation  along  with  accurate  position  measurements  from  onboard  GPS 
receivers,  accelerometer  measurements,  satellite  attitude  changes  and  inter-satellite 
range  changes.  All  of  these  measurements  allow  GRACE  to  measure  the  finer  effects 
of  the  gravity  field  that  cannot  be  seen  with  an  individual  vehicle. 

Both  satellites  were  launched  on  17  March  2002  into  co-planar,  near  circular, 
polar  orbits  at  500  km  altitude;  the  separation  of  the  satellites  is  kept  between  170 
and  220  km,  which  requires  a  maneuver  approximately  every  30-60  days.  Other  than 
maintaining  the  separation  of  the  satellites,  the  satellites  are  left  to  the  natural  effects 
of  the  gravity  field  and  the  orbits  are  allowed  to  decay;  as  of  Aug  2007,  the  satellites 
semimajor  axes  had  decayed  approximately  27  km  [1], 

The  GRACE-based  Earth  Gravity  Model  (EGM)  is  published  approximately 
every  30  days.  The  accuracy  of  the  model  comes  from  the  individual  measurements 
of  each  satellite  as  well  as  the  comparison  of  the  measurements.  The  comparison  of 
the  satellite  measurements  is  how  the  finer  effects  are  determined;  since  the  satellites 
pass  over  nearly  the  same  position  on  Earth  at  slightly  different  times,  the  smallest 
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temporal  and  spatial  variations  of  the  gravity  field  will  show  np  in  the  differences  in 
the  satellite  measurements. 

The  measurements  of  the  GRACE  satellites  have  shown  significant  improvement 
in  the  EGM,  for  terms  below  about  110,  over  previous  models  such  as  EGM-96  [16]. 
The  GRACE  based  EGM  has  also  been  combined  with  EGM-96  to  create  a  full 
degree/order  360  model  (EIGEN-CG01C)  [16].  And  yet,  work  continues  to  bring 
further  improvement  to  the  models. 

2.2  Jason-1 

The  Jason-1  satellite  carries  on  the  mission  of  the  Topex/Poseiden  satellite. 
Radar  altimetry  is  used  to  measure  surface  height  of  the  world’s  ocean  to  an  accu¬ 
racy  of  3.3  cm  [3];  some  of  the  other  mission  objectives  include  understanding  ocean 
circulation  and  understanding  how  changes  in  the  oceans  and  the  atmosphere  are 
related. 

Jason-1  was  launched  on  7  December  2001  into  a  circular  orbit  at  an  inclination 
of  66°  and  a  mean  altitude  of  1336  km.  This  allows  for  visibility  of  90%  of  the  world’s 
oceans  with  a  nominal  repeat  period  of  10  days.  Maneuvers  are  scheduled  between 
repeat  periods  whenever  possible  to  ensure  the  most  continuity  in  measurements. 

A  key  to  obtaining  accurate  measurements  of  the  ocean  heights  is  knowing  the 
orbital  position  of  the  satellite  to  high  precision.  For  Jason-1,  the  orbital  positions  are 
obtained  from  a  combination  of  Satellite  Laser  Ranging  (SLR)  and  the  Doppler  Or¬ 
bitography  and  Radio-positioning  Integrated  by  Satellite  (DORIS)  tracking  systems. 
This  method  is  able  to  provide  position  data  that  is  accurate  to  within  3  cm  [3]. 

2.3  Modern  Orbit  Determination 

The  position  data  for  the  GRACE  and  Jason-1  satellites  are  quite  precise,  both 
being  within  a  few  centimeters.  Unfortunately,  these  are  only  obtained  after  the 
satellite  has  passed  through  the  given  positions.  Trying  to  predict  a  satellite’s  future 
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position  to  the  same  sort  of  accuracy  is  not  currently  possible  for  any  extended  length 
of  time. 

The  methods  that  are  used  today  to  do  orbit  determination,  prediction  and 
modeling  rely  on  numerical  integration  of  the  perturbed  2BP,  or  analytical  solutions  to 
the  equations  of  motion.  Numerical  integration  can  be  very  computationally  expensive 
and  has  a  limit  to  the  period  of  validity.  In  the  past  it  could  take  as  much  time  to 
numerically  integrate  an  orbit  as  it  took  for  the  satellite  to  move  through  the  orbit  to 
the  point  to  which  you  were  integrating.  As  the  speed  of  computers  has  increased,  the 
time  required  to  perform  these  integrations  has  decreased.  However,  at  the  same  time, 
the  number  of  orbits  needing  to  be  integrated  has  increased  and  the  computational 
load  remains  high. 

General  analytical  solutions  of  the  equations  of  motion  through  series  expan¬ 
sions  are  computationally  less  expensive  than  numerical  integration  of  the  orbit;  the 
solution  provides  the  osculating  orbital  elements  as  a  function  of  time,  which  are  used 
to  determine  future  positions  based  on  initial  conditions  [18].  Unfortunately,  there 
is  generally  a  high  level  of  complexity  in  the  series  expansions.  To  reduce  the  com¬ 
plexity  of  the  solution,  the  series  are  usually  truncated  after  the  first  few  terms.  This 
also  solves  the  problem  of  small  divisors,  which  is  common  when  dealing  with  orbits 
that  have  small  eccentricities  [17,18];  however,  it  introduces  small  error  sources  which 
build  up  in  the  solution  over  time. 

Current  orbit  determination  methods  are  also  restricted  by  how  far  into  the 
future  they  can  model  the  orbit  due  to  the  buildup  of  errors  in  the  accuracy.  Because  of 
the  high  computational  load  required  for  the  integrations,  the  number  of  perturbations 
included  is  kept  fairly  low.  For  example  the  Earth’s  gravitational  potential  has  been 
measured  to  degree  and  order  360,  however,  to  reduce  the  computational  load,  the 
J2  term  (which  is  the  largest  term  in  the  geopotential)  is  almost  always  included, 
but  depending  on  the  computational  load  that  can  be  handled,  there  may  be  no 
other  terms  included,  or  there  may  be  quite  a  few  other  terms  included.  In  the 
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analytical  solution,  this  problem  arises  from  the  truncation  of  the  series  expansions. 
The  magnitude  of  the  terms  beyond  J2  are  on  the  order  of  10-6  and  smaller,  so 
this  may  not  have  a  significant  effect  as  long  as  the  time  period  for  the  prediction 
is  limited,  but  as  the  time  increases  the  uncertainty  of  the  solution  increases  and  its 
validity  decreases. 

One  area  that  is  hindered  by  the  current  methods  of  orbit  determination  is  the 
formation  flying  of  satellites.  Each  satellite  in  a  formation  requires  its  own  orbit 
model,  which  leads  to  its  own  level  of  uncertainty.  Assuming  the  formation  needs 
to  stay  within  certain  tolerances  concerning  separation  and  geometry,  the  problem 
becomes  very  difficult.  A  more  accurate  orbit  determination  method  with  a  longer 
period  of  validity  is  needed  to  make  formation  flight  of  satellites  more  feasible. 

2.4  Kolmogrov,  Arnold,  Moser  Theorem 

The  KAM  Theorem  is  the  result  developed  by  Kolmogrov  [9],  then  proved  by 
Arnold  [2]  and  Moser  [14],  to  describe  the  motion  of  a  nearly  integrable  Hamiltonian 
system  that  is  only  subject  to  small,  smooth  perturbations.  The  Hamiltonian  for  such 
a  system  is  given  by 


He{I,v)  =  h0{I)+eh{I,ip),  (2.1) 

where  ha  and  h  are  real-analytic  functions,  e  is  a  small  real  parameter,  and  /  and  ip  are 
symplectic  action-angle  variables.  The  theory  states  that  such  a  system  will  lie  on  a 
torus  that  can  be  described  by  the  fundamental  frequencies  of  the  orbit  and  associated 
Fourier  series  coefficients;  an  N-dimensional  Hamiltonian  will  have  N  frequencies  and 
will  lie  on  a  torus  in  2N-dimensional  phase  space.  This  equality  of  the  number  of 
coordinates  and  frequencies  is  mandated  by  the  Hamiltonian- Jacobi  theorem. 

KAM  theorem  has  been  shown  to  approximate  well  the  Restricted,  Circular, 
Planar,  Three-Body  Problem  (RCPTBP),  specifically  the  motion  of  asteroids  under 
the  effects  of  the  Sun  and  Jupiter  [5-7].  Celletti  and  Chierchia  [6]  showed  that  the 
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asteroid  12-Victoria,  in  a  Keplerian  orbit  about  the  Sun  and  perturbed  by  Jupiter, 
could  be  modeled  using  the  KAM  theorem;  in  this  case  the  perturbation  is  equal  to 
the  Jupiter-Sun  mass  ratio,  £  ~  10~3.  They  were  able  to  show  that  the  results  of 
their  analysis  agreed  well  with  observed  data  for  the  system. 

More  recently,  Wiesel  [19]  has  shown  that  frequencies  for  Earth  orbits  could  be 
identified,  and  the  orbits  could  be  fit  to  KAM  Tori  [13,20].  Assuming  the  torus  is 
exact,  any  trajectory  that  lies  on  the  torus  will  remain  on  the  torus  in  the  future. 
This  differs  from  perturbation  theory,  which  is  only  an  approximation  of  what  will 
happen  in  the  future  and,  therefore,  must  be  updated  due  to  error  growth. 

2.5  Orbital  Dynamics 

To  apply  the  KAM  theorem  to  Earth  orbiting  satellites,  a  reference  frame  must 
be  chosen  where  the  torus  will  be  a  static  geometric  structure;  this  means  we  need  a 
coordinate  frame  that  rotates  with  the  Earth  [19].  The  Earth  Centered  Earth  Fixed 
(ECEF)  coordinate  frame  is  chosen  as  the  frame  of  reference  for  this  work,  and  is 
defined  with  the  x  component  through  the  Prime  Meridian  at  the  equator,  the  £ 
component  northward  through  the  axis  of  rotation  and  the  y  component  pointing  in 
the  right  hand  sense. 

In  the  ECEF  frame,  the  coordinates  of  the  Hamiltonian  are  the  positions  x, 
y,  and  z  and  the  momenta  are  equal  to  the  inertial  velocities  resolved  along  the 
coordinate  axes  [19].  The  momenta  are  given  in  Equation  2.2 


px  =  x-  ujey 
py  =  ij  +  u ;ex 

Pz  =  z  (2.2) 


The  Hamiltonian  is  found  from  H  =  YliPiQi  —  T  +  V]  rearranging  Equation 
(2.2),  to  get  the  time  derivatives  of  the  coordinates,  and  substituting  back  in  gives 


the  Hamiltonian  in  terms  of  the  coordinates  and  momenta  only.  The  resulting  Hamil¬ 
tonian  for  Earth  orbiting  satellites  is  given  by, 


x  (Cnm  cos  m\  +  Snm  sin  mX) 


(2.3) 


where  P®  is  the  equatorial  radius  of  Earth,  p,  is  the  gravitational  parameter,  cu®  is  the 
rotation  rate  of  the  Earth  and  Cnm  and  Snm  are  the  gravity  field  coefficients.  The  P™ 
are  the  associated  Legendre  Polynomials,  and  r,  S  and  A  are  the  radius,  geocentric 
latitude  and  east  longitude  of  the  of  the  satellite,  respectively,  and  are  defined  in 
equation  (2.4). 


sin  5  = 


\Jx2  +  y 2 


tan  A  =  — 

x 


(2.4) 


Notice,  in  this  reference  frame  the  Hamiltonian  has  no  dependence  on  time  and  is 
thus  a  constant  of  the  motion. 

2. 6  The  Geopotential 

For  Earth  orbiting  satellites,  the  perturbation  that  is  of  interest  in  the  applica¬ 
tion  of  the  KAM  theorem  is  the  geopotential,  which  is  due  to  the  non- homogeneity 
of  the  planet.  Figure  2.1  displays  the  differences  between  the  theoretical  geopotential 
of  a  smooth  spherical  Earth  and  the  actual  geopotential  as  measured  by  GRACE  [1], 
The  areas  where  the  potential  is  highest  when  compared  to  the  smooth  Earth  are 
raised  and  colored  red,  while  the  areas  of  lowest  potential  are  recessed  and  colored 
blue. 
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Figure  2.1:  GRACE  geopotential  differences  from  smooth  spher¬ 
ical  Earth  [1], 
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In  addition  to  the  geopotential,  Earth  orbiting  satellites  also  experience  atmo¬ 
spheric  drag,  solar  radiation  pressure  and  other  effects  depending  on  their  orbital 
altitude.  For  low  altitude  satellites  (below  ~  350  —  400  km)  atmospheric  drag  is  the 
dominant  perturbation.  As  the  orbital  altitude  nears  Geosynchronous  Orbit  (GEO) 
the  effects  of  the  gravity  held  drop  off  and  effects  from  solar  radiation  pressure  in¬ 
crease.  Between  these  two  regions  is  where  the  geopotential  dominates  and  where  this 
author  seeks  to  apply  the  KAM  theorem. 

Using  the  KAM  theorem  to  model  satellite  motion,  based  on  real  orbital  data, 
takes  into  account  all  of  the  geopotential  effects  on  the  satellite.  This  is  a  huge 
advantage  over  current  methods  of  orbit  determination  and  prediction  which  typically 
use  only  a  subset  of  the  full  geopotential.  Since  all  of  the  geopotential  terms  are 
included  in  a  KAM  torus  defined  by  the  analysis  of  on  orbit  data,  one  should  be  able 
to  predict  more  accurately  and  for  longer  periods  of  time  than  with  current  methods. 

Including  the  whole  of  the  geopotential  in  the  orbit  solution  may  also  allow  for 
modeling  of  other  perturbations.  Small  changes  in  drag  on  the  satellite,  effects  from 
the  Moon,  Sun  and  other  celestial  bodies  may  be  masked  or  partially  masked  by  the 
effects  of  the  geopotential,  especially  when  the  effects  are  on  the  order  of  the  smallest 
terms  in  the  geopotential.  Once  the  geopotential  is  removed  it  may  be  possible  to 
more  accurately  study  and  model  such  effects. 

2. 7  Summary 

Modern  methods  of  orbit  determination  and  prediction  are  hindered  by  not 
including  all  of  the  terms  of  the  geopotential  when  integrating  the  orbit,  which  leads 
to  errors  in  the  accuracy  that  build  as  the  period  of  the  integration  is  increased. 
The  KAM  theorem  may  provide  a  better  method,  using  actual  position  data  for  the 
analysis  to  End  the  fundamental  frequencies  and  the  corresponding  coefficients.  This 
work  attempts  to  apply  KAM  theorem  to  the  GRACE  and  Jason-1  orbits,  to  show 
that  the  orbits  could  be  fit  to  tori  for  the  purpose  of  modeling  and  prediction. 
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III.  Method 


This  chapter  begins  with  a  discussion  of  the  satellites  and  corresponding  data.  The 
fundamental  frequencies  are  defined,  followed  by  a  discussion  of  the  algorithm  used  to 
find  the  frequencies  from  the  data.  Finally,  the  construction  of  the  torus  is  explained, 
along  with  how  the  torus  is  compared  to  the  original  data  to  understand  the  accuracy 
of  defining  the  orbit  with  this  method. 

3.1  Data 

The  GRACE  data  hies  were  obtained  from  the  Physical  Oceanography  Dis¬ 
tributed  Active  Archive  Center  (PO.DAAC)  at  NASA’s  Jet  Propulsion  Laboratory. 
The  positions  are  obtained  from  the  onboard  GPS  receivers;  they  are  given  in  the 
ECEF  coordinate  frame  at  60  second  intervals  and  are  accurate  to  within  a  few  cen¬ 
timeters  [16].  The  data  hies  also  contain  position  errors,  velocities,  velocity  errors, 
data  validation  fields  and  other  fields  that  are  not  required  for  this  analysis  [4],  A 
Matlab®  script  was  written  to  extract  the  x,  y,  and  z  positions  from  the  data  hies. 

The  Jason-1  data  hies  were  obtained  from  NASA’s  Goddard  Space  Flight  Center 
(GSFC).  The  positions  are  determined  from  a  combination  of  the  SLR  and  DORIS 
tracking  of  the  satellite;  they  are  given  in  the  ECEF  coordinate  frame  at  60  second 
intervals  and  are  accurate  to  within  3-4  cm  [3].  The  data  hies  contain  10-day  intervals 
that  correspond  to  the  10-day  repeat  cycle  of  the  orbit.  To  ensure  enough  data  was 
available  for  the  analysis,  three  repeat  cycles  were  used.  Cycles  133-135  were  chosen 
because  they  followed  a  maneuver;  this  was  done  to  reduce  the  likelihood  that  a 
maneuver  would  occur  during  the  period  of  analysis.  Again,  the  data  hies  contained 
more  information  than  was  needed,  so  a  Matlab®  script  was  written  to  extract  the 
x,  y,  and  z  positions  from  the  data  hies. 

The  position  data  sets  for  both  GRACE  and  Jason- 1  are  on  the  order  of  ~ 
100,000  meters  and  are  separated  by  60  seconds.  This  author  uses  the  dimension¬ 
less  quantities  of  Distance  Units  (DU)  and  Time  Units  (TU)  during  analysis,  which 
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reduces  the  data  to  on  the  order  of  ~  1.2  DU. 


1  DU  =  6378136.3  m 


1  TU  =  806.810988  s 


(3.1) 

(3.2) 


In  addition  to  reducing  the  magnitudes  of  the  position  data,  the  choice  of  DUs  and 
TUs  results  in  the  radius  of  Earth,  R®,  and  the  gravitational  parameter,  /z,  being 
equal  to  1.  This  simplifies  the  calculation  of  the  estimate  frequencies  in  Equations 

(3.3)-(3.5)  below. 

3.2  Fundamental  Frequencies 

An  Earth  orbiting  satellite  will  have  three  fundamental  frequencies:  the  anoma¬ 
listic  frequency,  the  Earth’s  rotation  rate  plus  the  nodal  regression  rate,  and  the 
apsidal  regression  rate.  These  frequencies  determine  the  time  is  takes  for  the  satellite 
orbit  to  traverse  the  dimensions  of  the  torus.  These  frequencies  are  approximated  by 
Equations  (3.3)-(3.5),  which  include  only  the  J2  effects  on  the  mean  motion,  the  right 
ascension  of  the  ascending  node,  and  the  argument  of  perigee  [17]. 


(3.3) 


(3.4) 


(3.5) 
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The  total  time  it  takes  for  a  satellite  to  traverse  the  entire  torus  is  dependent  on  the 
smallest  frequency.  It  gives  the  toroidal  period  as 

2n 

Tram  =  — •  (3.6) 

U>3 


3.3  Frequency  Determination  Algorithm 

Since  the  data  used  in  this  analysis  covers  a  finite  period  of  time,  it  seems 
natural  to  use  a  finite  time  Fourier  transform  method  of  the  form  in  Equation  3.7.  In 
this  work,  the  function  fit)  is  the  position  coordinates  for  the  orbit. 


1  f  me^xm  -  1  J  f(t)e-"‘X(t)dt 


(3.7) 


The  Fourier  transform  algorithm  used  in  this  work  is  based  on  Laskar’s  [10-12] 
Numerical  Algorithm  of  the  Fundamental  Frequency  (NAFF),  which  has  been  shown 
to  identify  the  frequencies  of  a  quasiperiodic  function  to  a  higher  precision  than  a 
simple  Fast  Fourier  Transform  (FFT).  An  FFT  typically  assumes  the  function  has  a 
period  of  2 r  over  a  range  of  [— r  :  r],  with  an  accuracy  proportional  to  1/r;  Laskar’s 
NAFF  does  not  make  this  assumption  and,  with  the  use  of  a  window  function,  y, 
converges  on  the  frequencies  with  an  accuracy  proportional  to  1/r4.  The  NAFF  finds 
the  frequency  by  iterating  over  Equation  (3.8). 

=  \  J  f  / T)dt  (3.8) 

The  Hanning  window  function  is  used  in  this  analysis  and  is  given  by  Equation  (3.9). 

y(t/r)  =  1  +  cos  (pp\  (3.9) 
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The  approximate  frequencies  from  Equations  (3.3)  -  (3.5)  are  input  into  the  algorithm, 
as  initial  guesses,  to  reduce  the  number  of  iterations.  As  noted  above,  the  precision 
of  the  frequency  found  is  higher  than  it  would  be  with  a  simple  FFT.  However,  this 
doesn’t  mean  that  the  frequencies  are  exact,  just  very  close. 

It  should  also  be  noted  that  the  Hanning  window  function  accelerates  the  con¬ 
vergence  of  finding  the  peak  frequency;  higher  order  Hanning  functions  result  in  faster 
convergence,  though  Laskar  shows  that  the  cost-to-benefit  drops  off  after  windows  of 
order  3-5  [11].  Higher  order  functions  also  result  in  wider  peaks,  which  leads  to 
more  uncertainty  in  the  frequency.  This  author  uses  the  2nd  order  Hanning  window 
function,  which  seems  to  work  fairly  well. 

Figure  3.1  shows  the  results  of  the  Fourier  transform  of  the  GRACE  data.  The 
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Figure  3.1:  PSD  plot  for  the  GRACE-A  satellite  showing  how 
the  modified  Laskar  algorithm  finds  the  frequencies. 
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triplet  structure  of  the  frequencies  is  due  to  the  fact  that  u>2,  the  combination  of  cu® 
and  the  nodal  regression  rate,  is  seen  in  the  x  and  y  components  of  the  position, 
but  it  does  not  affect  the  z.  This  can  be  explained  by  first,  noting  that  all  of  the 
components  change  with  the  frequency  of  the  orbit,  uq,  but  the  x  and  y  components 
also  change  with  cu®.  This  is  highlighted  by  the  peak  in  the  x  and  y  frequencies,  at  ~ 
0.0585  rad/TU,  which  corresponds  to  cu®,  but  no  such  peak  exists  in  the  z  frequency. 
Second,  the  regression  of  the  node  is  a  rotation  about  the  z  axis,  so  it  won’t  show  up  in 
the  2;  component  of  the  orbit.  The  triplet  structure  of  the  frequencies  is  also  repeated 
at  equal  intervals,  which  correspond  to  integer  multiples  of  the  orbital  frequency, 
nuii.  The  smaller  peaks  that  flank  each  large  peak  are  higher  order  combinations  of 
nuj\  ±  mw2,  where  n  and  m  denote  integer  multiples  of  the  frequencies  [15]. 

The  width  of  the  peaks  is  partially  a  consequence  of  the  Hanning  window  and 
partially  a  consequence  of  the  finite  span  of  data.  Using  a  larger  span  of  data  would 
reduce  the  width  of  the  peaks,  possibly  resulting  in  more  accurate  identification  of  the 
frequencies.  However,  the  issue  that  one  encounters  when  analyzing  longer  periods  of 
on  orbit  satellite  data,  is  maneuvering.  If  a  torus  were  already  created  for  an  orbit,  a 
maneuver  would  change  the  frequencies,  the  series  coefficients,  or  both,  which  would 
require  the  torus  to  be  recalculated.  A  maneuver  during  the  analysis  of  the  data  could 
also  change  the  frequencies,  resulting  in  the  apparent  identification  of  multiple  peaks 
for  each  frequency,  which  would  invalidate  the  analysis. 

3.4  The  KAM  Torus 

After  the  frequencies  are  identified,  the  torus  can  be  formed  from  the  funda¬ 
mental  frequencies  and  the  Fourier  coefficients.  The  Fourier  coefficients  are  defined 
for  a  single  frequency  line  in  [20]  as, 

C  =  23L4(w),  5  =  2AA(w)  (3.10) 
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where  A(cu)  is  the  peak  amplitude  of  the  frequency  multiples  found  during  the  Fourier 
transform  analysis.  With  the  frequencies  and  the  coefficients,  the  coordinates  of  the 
KAM  torus  can  be  calculated  by  Equation  (3.11);  j  is  the  series  limit  of  the  coefficients, 
j  is  a  summation  vector,  and  cu  is  the  vector  of  torus  basis  frequencies  [20] . 

q(t)  =  Co  +  cos  (j  •  cut)  +  Sj  sin(j  •  cut)}  (3.11) 

j 

In  this  work,  a  series  limit  of  10  is  used  when  calculating  the  Fourier  coefficients;  this 
means  the  coefficients  are  defined  for  the  frequency  harmonics  up  to  lOoqilOo^.  The 
constant  term  Cq  is  found  from  the  initial  conditions. 

Using  a  known  starting  position,  the  positions  for  any  point  in  the  future  (or 
past)  can  then  be  calculated  and  compared  to  the  real  data  to  get  the  residuals  and 
determine  how  well  the  torus  approximates  the  orbit.  The  torus  can  have  errors  in 
both  the  frequencies  and  the  Fourier  series  coefficients,  but  they  show  up  in  distinct 
and  different  ways.  If  there  are  small  errors  in  the  frequencies,  they  will  show  up  as 
linear  error  growth  over  time.  The  coefficients  define  the  shape  and  position  of  the 
torus  in  space,  so  errors  show  up  in  the  torus  as  either  a  shift  from  where  it  should 
be  in  space,  or  in  a  shape  difference  from  that  of  the  true  torus. 

3. 5  Summary 

The  fundamental  frequencies  for  the  GRACE  and  Jason-1  orbits  are  found  from 
the  modified  Laskar  algorithm  analysis  of  the  precise  orbital  position  data.  The 
Fourier  coefficients  are  then  calculated  to  complete  the  definition  of  the  KAM  tori  for 
each  satellite,  allowing  the  orbit  to  be  modeled.  The  positions  are  then  calculated 
from  the  torus  and  compared  to  the  real  position  data  to  determine  how  well  the 
torus  models  the  orbit. 
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IV.  Results  and  Findings 

This  chapter  shows  that  discrete  frequencies  exist  for  both  the  GRACE  and  Jason- 1 
satellites.  These  are  then  used  to  create  the  tori.  The  tori  are  used  to  predict  the 
positions  of  the  satellites  and  compared  to  the  real  orbital  data. 

4.1  GRACE 

This  section  will  cover  the  frequency  estimates  for  the  GRACE-A  satellite,  fol¬ 
lowed  by  the  frequencies  found  using  the  modified  Laskar  algorithm.  Finally,  the 
torus  is  compared  with  the  real  orbital  data. 

4.1.1  Frequency  Estimates.  The  estimates  for  the  GRACE  frequencies 
are  calculated  from  Equations  (3.3)-(3.5)  and  are  given  in  Equations  (4.1)-(4.3). 
The  values  used  for  calculating  the  frequencies  are  given  as  a  ~  6843000  m,  n  — 
3986004. 415E+8  m3/s2,  R#  =  6378136.3  m,  i  «  89°,  e  «  .0017,  J2  =  4.84165339915E- 
4,  and  u ;e  =  7.292115857916E  -  5  rad/s. 

(hi  =  0.899567967941131  rad/TU  (4.1) 

=  1.114967418186E  -  3  rad/s 

oj2  =  5.8843500509237E  -  2  rad/TU  (4.2) 

=  7.2933439658E  -  5  rad/s 

Cu3  =  2.83439973794 E  -  4  rad/TU  (4.3) 

=  3.51309014E  -  7  rad/s 

4-1.2  Laskar  Frequency  Output.  The  Laskar  algorithm  only  identified  the 
first  two  frequencies  for  the  GRACE  orbit;  a  value  for  the  third  frequency  was  given, 
but  it  was  on  the  order  of  10~30,  which  is  effectively  zero  and  assumed  to  be  erroneous. 
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This  is  not  altogether  surprising,  since  the  orbit  is  very  circular,  e  ~  .0017465.  This 
problem  occurs  for  two  reasons.  First,  based  on  the  approximate  frequency  from 
equation  (4.3),  above,  we  expect  the  period  of  a  torus  for  this  orbit  to  be  ~  207 
days.  Only  31  days  worth  of  data  were  used  (1-31  Mar  08),  so  the  perigee  only  moved 
.721  rad ;  this  makes  it  difficult  to  see  the  movement  of  the  perigee,  given  the 
amount  of  data  used.  Second,  as  mentioned  in  Section  3.3,  the  peak  that  identifies 
the  frequency  has  a  width.  Since  the  frequency  is  already  close  to  zero,  the  width 
of  the  peak  identified  for  the  smallest  frequency  may  cross  zero,  which  leads  to  the 
error. 

The  first  two  frequencies  are  found  to  be  very  close  to  the  approximate  frequen¬ 
cies  given  by  Equations  (4.1)  and  (4.2). 


uq  =  0.899330102615437  rad/TU 


(4.4) 


u2  =  5.890847214269741E  -  2  rad/TU 


(4.5) 


Figures  4. 1-4.4  are  Power  Spectral  Density  (PSD)  plots  for  the  GRACE-A  data. 
The  largest  peak  in  Figure  4.1  corresponds  to  aq  and  appears  in  the  analysis  of  the  z 
component  because  it  is  unaffected  by  the  Earth’s  rotation  and  nodal  regression  rates. 
Alternatively,  x  and  y  are  affected  by  the  Earth’s  rotation  and  the  nodal  regression 
and  so,  the  largest  peaks  in  Figure  4.2  and  Figure  4.3  come  from  uq  ±  uq. 

Figure  4.4  has  all  three  components  plotted  to  show  how  the  peaks  in  the  x 
and  y  components  straddle  the  z  peak.  The  smaller  frequency  range  displayed,  also 
shows  the  width  of  the  peaks  more  clearly  than  Figures  4.1  -  4.3.  As  discussed  in 
Section  3.3,  the  width  of  the  peaks  is  a  source  of  error  in  the  frequencies.  This  may 
in  fact  be  the  reason  uq  was  not  identified.  In  the  z  frequency  there  is  a  spike  right  at 
zero.  It  is  possible  this  spike  corresponds  to  uq,  but  due  to  the  width  the  top  of  the 
peak  cannot  be  identified,  which  makes  it  impossible  to  find  the  frequency.  Further 
analysis  is  required  to  confirm  this  hypothesis. 
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Figure  4.1:  z  frequency  PSD  plot  for  the  GRACE-A  satellite. 
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Figure  4.2:  x  frequency  PSD  plot  for  the  GRACE-A  satellite. 
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Figure  4.3:  y  frequency  PSD  plot  for  the  GRACE-A  satellite. 
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Figure  4.4:  GRACE-A  PSD  plot  highlighting  the  triplet  struc¬ 
ture  of  the  frequencies. 
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4-1.3  Comparing  KAM  to  Real  Data.  After  the  frequencies  were  identified 
by  the  Laskar  algorithm,  the  torus  is  formed  and  compared  to  the  real  data,  and  the 
residuals  were  calculated.  The  initial  residuals,  which  are  displayed  in  Figure  4.5,  are 
nearly  zero  at  the  initial  conditions,  t  —  0.  However,  frequency  errors  are  present  in 
both  frequencies  as  evidenced  by  the  error  growth  as  t  goes  toward  ±tfinai. 
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Figure  4.5:  This  shows  the  residual  growth  over  time  using  the 
frequencies  found  from  the  Laskar  algorithm. 


’xresidual.txt’ 

’yresidual.txt’ 


The  frequency  errors  must  be  small  since  they  show  up  as  linear  error  growth; 
the  explanation  is  given  by 

X  =  C  cos(edi  +  5u±)t  +  S  sin(edi  +  Sui)t 
=  C (cos  oj\t  cos  5uj\t  —  sin  e opt  sin  Suit) 

+  ^(sin  uit  cos  Suit  +  cos  Uit  sin  Suit) 

~  C  cos  edit  —  C Suit  sin  edit  +  S  sin  edit  +  SSuitcosuit  (4.6) 
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using  standard  trigonometric  identities  and  the  small  angle  approximation.  The  <5uqt 
terms  that  are  left,  result  in  the  linear  error  growth  visible  in  Figure  4.5;  the  same 
result  can  be  shown  for  5u2- 

From  the  slope  of  the  residual  peaks,  we  can  approximate  the  magnitude  of 
the  frequency  errors.  Since  the  torus  only  requires  the  frequencies  and  the  Fourier 
coefficients,  it  is  possible  to  adjust  the  frequencies,  without  having  to  rerun  the  entire 
frequency  analysis,  to  try  and  reduce  the  error.  While  adjusting  the  frequencies,  the 
Fourier  coefficients  are  assumed  to  be  correct  and  held  constant. 

The  frequency  errors  for  uq  and  uq  were  found  to  be  on  the  order  of  10~5  rad/TU. 
Adjusting  only  uq,  the  frequency  error  was  found  to  be  ~  7 77  —  5  rad/TU.  Then, 
adjusting  only  uq,  the  frequency  error  was  found  to  be  ~  5 E  —  5  rad/TU ;  further 
adjustments  were  made  by  making  smaller  changes  to  the  frequencies.  Figure  4.6 
shows  how  the  residuals  improved  after  the  frequencies  were  adjusted. 


Figure  4.6:  Residual  improvements  with  5uq  =  7.477—5  rad/TU 
and  <5uq  =  5.477  —  5  rad/TU. 
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Further  adjustments  revealed  a  limit  to  the  improvement  that  could  be  gained 
by  the  frequency  changes.  Figure  4.7  shows  the  residuals  for  the  frequencies  adjust 
by  Suji  =  7.7 E  —  5  rad/TU  and  Su>2  =  5.7 E  —  5  rad/TU .  Over  all,  the  residuals  are 
much  improved  over  the  initial  findings.  However,  the  residuals  have  begun  to  show 
a  quadratic  form  in  the  £  component.  This  is  likely  due  to  atmospheric  drag,  which 
appears  as  quadratic  changes  in  the  mean  anomaly  of  the  orbit  [17].  Atmospheric  drag 
is  not  modeled  by  the  KAM  theorem,  and  must  be  handled  by  some  other  means. 
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Figure  4.7:  Residual  improvements  with  Su q  =  7.7E—5  rad/TU 
and  5u>2  =  5.7 E  —  5  rad/TU. 
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Looking  closer  at  the  area  where  the  residuals  are  the  smallest,  Figure  4.8  shows 
they  are  still  fairly  high.  The  z  coordinate  residuals  are  ~  20  km ,  and  the  x  and  y 
coordinate  residuals  are  off  by  as  much  as  40  —  90  km. 

Some  of  the  error  still  present  in  the  residuals  is  due  to  errors  in  the  coefficients. 
This  is  most  likely  the  case  with  the  offset  that  is  present  in  each  of  the  components; 
the  y  component  shows  the  largest  offset,  but  all  three  have  some  level  of  offset.  Figure 
4.8  also  shows  a  distinct,  repeating  pattern  in  the  x  and  y  residuals.  Looking  more 
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Figure  4.8:  This  shows  the  error  still  present  in  the  region  where 
the  frequency  adjustment  is  most  effective. 


closely  at  the  x  residuals,  Figure  4.9  shows  that  the  frequency  of  the  peaks  is  nearly 
equal  to  the  orbital  frequency.  There  are  also  consistently  8  peaks  in  each  repetition. 
This  suggests  there  may  be  a  dominate  error  in  the  coefficient  defined  by  the  eighth 
harmonic  of  uj\ .  However,  further  analysis  is  needed  to  confirm  this  and  resolve  the 
issue. 

Unfortunately,  these  remaining  errors  are  too  large  for  the  torus  to  be  useful  for 
modeling  the  orbit.  If  these  errors  could  be  driven  down  further,  through  refinement  of 
the  torus,  it  may  be  possible  to  model  the  GRACE  orbit  with  a  KAM  torus.  However, 
perturbation  theory  would  still  be  required  to  account  for  the  effects  of  atmospheric 
drag. 
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4-2  Jason-1 

This  section  will  cover  the  frequency  estimates  for  the  Jason-1  satellite,  followed 
by  the  frequencies  found  using  the  modified  Laskar  algorithm.  Finally,  the  torus  will 
be  compared  with  the  real  orbital  data. 

4-2.1  Frequency  Estimates.  The  estimates  for  the  Jason-1  frequencies 
are  calculated  from  Equations  (3.3)-(3.5)  and  are  given  by  Equations  (4.7)-(4.9). 
The  values  used  for  calculating  the  frequencies  are  given  as  a  =  7715636.3  m, 
fi  =  3986004.415E  +  8  m3/s2,  Re  =  6378136.3  m,  i  «  66°,  e  =  .0007465,  J2  = 
4.8416533991517  -  4,  and  a;©  =  7.292115857916E  -  5  rad/s. 

u>i  =  0.751499356368578  rad/TU  (4.7) 

=  9.3156093090217  -  4  rad/s 
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cJ2  =  5.8985058639173E  -  2  rad/TU 
=  7.3108893552 E  -  5  rad/s 


(4.8) 


a)3  =  3.273524554577  -  5  rad/TU  (4.9) 

=  4.057362477  -  8  rad/s 

4-2.2  Laskar  Frequency  Output.  The  Laskar  algorithm  only  identified  the 
first  two  frequencies  for  the  Jason-1  orbit;  a  value  for  the  third  frequency  was  given, 
but  it  was  on  the  order  of  10~28,  which  is  effectively  zero  and  considered  erroneous. 
This  is  not  altogether  surprising,  since  the  orbit  is  very  circular,  e  =  .0007465.  This 
problem  occurs  for  two  reasons.  First,  based  on  the  approximate  frequency  from 
equation  (4.9),  above,  we  expect  the  period  of  a  torus  for  this  orbit  to  be  ~  5  years. 
Only  as  30  days  worth  of  data  were  used,  so  the  perigee  only  moved  ~  .103  rad ; 
this  means  it  would  be  difficult  to  see  the  movement  of  the  perigee  in  the  period 
analyzed.  Second,  as  mentioned  in  section  3.3,  the  peak  that  identifies  the  frequency 
has  a  width.  Since  the  frequency  is  already  close  to  zero,  the  width  of  the  peak  for 
the  smallest  frequency  may  cross  zero,  which  leads  to  the  error. 

The  first  two  frequencies  are  found  to  be  very  close  to  the  approximate  frequen¬ 
cies  given  by  equations  (4.7)  and  (4.8). 

oq  =  0.751485615955426  rad/TU  (4.10) 

u2  =  5.91535231579612177  -  2  rad/TU  (4.11) 

Figures  4.10-4.13  are  PSD  plots  for  the  Jason-1  data.  The  largest  peak  in 
Figure  4.10  corresponds  to  uq.  As  discussed  above,  it  appears  in  the  analysis  of  the 
z  component  because  it  is  unaffected  by  the  Earth’s  rotation  and  nodal  regression 
rates.  Again,  x  and  y  are  affected  by  the  Earth’s  rotation  and  the  nodal  regression 
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Figure  4.10:  The  z  frequency  PSD  plot  for  the  Jason-1  satellite. 


and  the  largest  peaks  in  Figure  4.11  and  Figure  4.12  come  from  uq  iuq.  Figure  4.13 
has  all  three  components  plotted  to  show  how  the  peaks  in  the  x  and  y  components 
straddle  the  z  peak. 

Another  result  to  note  in  Figures  4.10-4.12  is  the  magnitudes  of  the  peaks.  In 
Fourier  analysis  of  quasiperiodic  systems,  one  expects  that  the  power  will  fall  off  with 
increasing  integer  multiples  of  the  frequencies  [15],  e.g.  A(aq)  >  A(2ui)  >  A(3uq). 
This  is  not  the  case  for  Jason-1.  Looking  specifically  at  Figure  4.10,  it  can  be  seen 
that  the  power  of  the  peak  for  3uq  is  less  than  the  peak  for  cq,  as  expected,  but 
larger  than  the  power  of  the  peak  for  2uq,  which  is  contrary  to  what  we  expect.  This 
suggests  there  may  be  an  error  in  the  coefficient  calculated  for  the  2uq  frequency  line, 
which  will  show  up  as  errors  in  the  residuals  when  the  torus  is  compared  to  the  real 
data. 
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Figure  4.11:  The  x  frequency  PSD  plot  for  the  Jason-1  satellite. 
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Figure  4.13:  Jason-1  PSD  plot  showing  the  triplet  structure  of 
the  frequencies. 

4-2.3  Comparing  KAM  to  Real  Data.  After  the  frequencies  were  identified 
by  the  Laskar  algorithm,  the  torus  is  formed  and  compared  to  the  real  data,  and  the 
residuals  were  calculated.  The  initial  residuals,  which  are  displayed  in  Figure  4.14, 
are  nearly  zero  at  the  initial  conditions,  t  =  0.  However,  a  frequency  error  in  ui2 
appears  as  error  growth  when  t  goes  toward  ±tfinaf,  the  lack  of  frequency  error  in  uq 
is  shown  by  the  fact  that  the  z  residuals  remain  near  zero  over  the  entire  time  span. 

The  frequency  error  is  small,  resulting  in  linear  error  growth  (see  equation  (4.6)). 
From  the  slope  of  the  residual  peaks,  we  can  approximate  the  magnitude  of  the  fre¬ 
quency  error.  Since  the  torus  only  requires  the  frequencies  and  the  Fourier  coefficients, 
it  is  possible  to  adjust  the  frequency,  without  having  to  rerun  the  entire  frequency 
analysis,  to  try  and  reduce  the  error.  While  adjusting  the  frequency  that  was  in  error, 
the  other  frequency  and  the  Fourier  coefficients  are  assumed  to  be  correct  and  held 
constant. 
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Figure  4.14:  This  shows  the  residual  growth  over  time  using  the 
frequencies  found  from  the  Laskar  algorithm. 

For  Jason-1,  the  frequency  error  was  found  to  be  on  the  order  of  10~5  rad/TU. 
After  some  analysis,  it  was  determined  that  the  error  was  ~  1.86177  —  5  rad/TU ,  with 
the  adjusted  frequency  given  in  Equation  (4.12). 

u2  =  5.917213315796121E  -  2  rad/TU  (4.12) 

Figure  4.15  shows  the  residuals  after  the  frequency  correction.  The  root  mean  square 
of  the  residuals  is  less  than  240  meters,  and  at  no  point  in  the  period  of  analysis  do 
the  residuals  exceed  1  km. 

Based  on  the  discussion  in  Section  4.2.2,  a  large  part  of  the  remaining  error 
is  believed  to  be  due  to  coefficient  errors.  Figure  4.16  highlights  the  x  residuals  for 
Jason-1,  which  do  not  show  a  periodic  repetition  over  the  span  of  data  like  that  seen 
in  the  GRACE  data.  This  suggests  that  unlike  the  GRACE  results,  where  it  seemed 
the  majority  of  the  coefficient  error  was  due  to  one  harmonic,  the  Jason-1  residual 
errors  may  be  due  to  multiple  coefficient  errors. 
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Figure  4.15:  This  shows  the  residuals  over  time  after  the  fre¬ 
quency  error  was  corrected. 
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4-3  Summary  of  Results 

KAM  tori  were  created  for  both  satellites,  but  after  correcting  for  frequency 
errors,  only  the  Jason-1  torus  was  able  to  produce  residuals  with  a  reasonable  level 
of  accuracy.  The  GRACE  data  produced  results  indicating  that  atmospheric  drag 
cannot  be  ignored,  which  hinders  the  ability  of  KAM  theorem  to  define  the  orbit 
without  further  studies.  The  Jason-1  torus  is  promising,  but  further  refinement  is 
needed  to  justify  its  use  for  orbit  modeling. 
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V.  Conclusions 


This  chapter  presents  the  conclusions  that  are  based  on  the  results  from  Section  IV. 
As  with  the  results,  each  satellite  is  discussed  individually.  Finally,  recommendations 
are  made  for  further  study  into  the  applicability  of  using  KAM  theorem  to  define  the 
orbits  for  GRACE  and  Jason-1. 

5.1  GRACE 

The  results  from  Section  4.1  suggest  that  the  GRACE  satellites  may  lie  on  a 
torus;  however,  the  effects  of  atmospheric  drag  cannot  be  ignored.  There  also  appear 
to  be  errors  in  the  coefficients  that  must  be  resolved  in  order  to  accurately  define 
the  torus.  If  the  coefficient  errors  can  be  resolved,  it  should  be  possible  to  define  the 
torus  and  use  perturbation  theory  to  account  for  the  effects  of  the  atmospheric  drag; 
however,  further  study  is  required  to  confirm  this. 

5.2  Jason-1 

The  results  from  Section  4.2  indicate  that  the  Jason-1  satellite  does  lie  on  a 
torus.  The  comparison  of  the  torus  data  with  the  real  data  shows  good  agreement, 
but  there  is  still  error  that  needs  to  be  removed.  Analysis  of  the  Fourier  coefficients 
should  be  able  to  resolve  the  errors  present. 

5.3  Recommendations  for  Further  Study 

The  methods  employed  in  this  work  provide  a  good  starting  point  for  defining 
KAM  tori  for  on  orbit  satellites.  The  next  step  should  take  the  frequencies  and  Fourier 
coefficients  identified  by  this  method  and  refine  them  further,  with  a  fitting  algorithm, 
such  as  least  squares. 

The  period  of  validity  also  needs  to  be  studied.  Characterization  of  how  the 
errors  continue  to  change  and  grow  outside  of  the  period  of  data  used  in  the  analysis 
is  needed.  The  Jason-1  torus  defined  in  this  study  could  be  used  for  a  first  analysis  of 
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this  problem,  but  a  torus  that  is  further  refined,  such  as  with  the  methods  mentioned 
above,  would  be  more  advantageous. 

Characterization  of  the  altitude  range  for  the  KAM  tori  should  also  be  done  to 
understand  where  the  non-geopotential  perturbations  can  be  ignored.  As  shown  by 
the  GRACE  results,  lower  altitude  satellites  must  deal  with  the  effects  of  atmospheric 
drag;  similarly,  much  higher  altitude  satellites  are  significantly  affected  by  the  Moon, 
Sun  and  solar  radiation  pressure,  and  must  take  into  account  these  effects.  Under¬ 
standing  where  these  regions  reside  will  allow  for  studies  into  how  the  perturbations 
affect  the  KAM  tori  and  how  the  effects  can  be  overcome. 
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Appendix  A.  Data  files 


A.l  GRACE 

The  GRACE  data  files  are  extracted  from  the  larger  Level-IB  (LIB)  data  hies 
that  are  available  from  the  PO.DAAC  at  JPL.  The  LIB  hies  are  produced  for  each 
day  and  can  produce  a  large  number  of  smaller  data  hies.  Table  A.l  describes  the 
format  of  the  LIB  data  hie  containing  the  positions  used  in  this  work. 

Table  A.l:  The  NASA  code  obtained  from  JPL  produces  the  GNV1B 
data  hies  containing  the  following  helds  [4]. 


Parameter 

Definition 

Data  Type 

Byte  Length 

Units 

gps.time 

seconds  past  noon  01- Jan- 2 000 

Integer 

4 

s 

GRACE id 

GRACE  satellite  identifier 

Character 

1 

N/A 

coord_ref 

Coordinate  reference  frame  where 

E  =  Earth  Fixed 

T  =  Inertial 

Character 

1 

N/A 

xpos 

Position,  x  value  (ITRF) 

Double 

Precision 

8 

m 

ypos 

Position,  y  value  (ITRF) 

Double 

Precision 

8 

m 

zpos 

Position,  z  value  (ITRF) 

Double 

Precision 

8 

m 

xpos_err 

Formal  error  on  x  position 

Double 

Precision 

8 

m 

ypos_err 

Formal  error  on  y  position 

Double 

Precision 

8 

m 

zpos_err 

Formal  error  on  z  position 

Double 

Precision 

8 

m 

xvel 

Velocity  along  x-axis  (ITRF) 

Double 

Precision 

8 

m 

yvel 

Velocity  along  y-axis  (ITRF) 

Double 

Precision 

8 

m 

zvel 

Velocity  along  z-axis  (ITRF) 

Double 

Precision 

8 

m 

xveLerr 

Formal  error  in  velocity  along  x-axis 

Double 

Precision 

8 

m 

yveLerr 

Formal  error  in  velocity  along  y-axis 

Double 

Precision 

8 

m 

zveLerr 

Formal  error  in  velocity  along  z-axis 

Double 

Precision 

8 

m 

qualflg 

Data  quality  flag  (LSB  =  bit  0) 

bit  0  =  Not  Defined 
bit  1  =  Not  Defined 

bit  2  =  overlap  data  missing  before  start  midnight 
bit  3  =  overlap  data  missing  after  start  midnight 
bit  4  =  overlap  data  missing  before  end  midnight 
bit  5  =  overlap  data  missing  after  end  midnight 
bit  6  =  Not  Defined 
bit  7  =  Not  Defined 

Unsigned 

Character 

1 

N/A 
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A.  2  Jason- 1 


The  Jason-1  data  files  are  available  from  the  GSFC  FTP  site  in  the  format 
described  in  Table  A. 2.  However,  to  ensure  the  data  is  accurate  it  must  be  run 
through  the  Hermite  interpolator,  also  available  from  the  GSFC  FTP  site  [8].  After 
the  interpolator  has  been  run  on  the  data,  it  outputs  the  positions  in  the  ECEF 
reference  frame,  as  well  as  velocities,  timing  data  and  other  information  that  was  not 
used  by  this  author. 

Table  A. 2:  The  Hermite  interpolator  code  obtained  from  NASA  ingests 
the  following  information,  to  produce  the  ECEF  position 
used  in  this  work. 


Record 

Item 

Format 

Description  of  Format 

N 

i 

D22.16 

Epoch  year,  month,  day,  hour,  min 
satellite  time  in  form  YYMMDDhhmm  (UTC) 

N 

2 

D22.16 

Epoch  seconds  satellite  time  (UTC) 

N 

3 

D22.16 

Sidereal  time  -  Greenwich  hour 
angle  satellite  epoch  (DEGREES) 

N 

4 

D22.16 

Polar  motion  X  (milli  arc-seconds) 

N 

5 

D22.16 

Polar  motion  Y  (milli  arc-seconds) 

N 

6 

D22.16 

Epoch  time  in  ET  days  from  Jan  0.0 
of  the  reference  year  of  the  arc  (days) 

N+l 

7-9 

3D22.16 

Satellite  Inertial  True  of  Date  (ITOD) 

X,Y,Z  elements  (meters) 

N+l 

10-12 

3D22.16 

Satellite  ITOD  X,Y,Z  elements  (meters/second) 

N+2 

13-15 

3D22.16 

Satellite  Earth  Centered  Fixed  (ECF) 

X,Y,Z  elements  (meters) 

N+2 

16-18 

3D22.16 

Satellite  ECF  X,Y,Z  elements  (meters/second) 

N+3 

19 

2211 

Flags  indicating  orbit  mode  (0=no,  l=yes) 

(1)  OCCULTATION  (0=sun,  l=shadow) 

Yaw  steering  events  (nominal  Beta’  example) 

(2)  +ON  (  15  to  80  DEG  BETA’) 

(3)  +OFF  (  0  to  15  DEG  BETA’) 

(4)  +HIGH  OFF  (  80  to  90  DEG  BETA’) 

(5)  +RAMP  UP  (  15  to  15.1DEG  BETA’) 

(6)  +RAMP  DOWN  (  15  to  14.9DEG  BETA’) 

(7)  +FLIP  (  +0  to  -0  DEG  BETA’) 
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Appendix  B.  Additional  GRACE  Results 

This  appendix  provides  results  for  the  GRACE  satellite  that  were  obtained  during 
the  analysis.  A  brute  force  method  for  changing  the  frequencies  is  highlighted  here, 
and  was  used  to  arrive  at  the  final  results  and  conclusions  presented  in  main  body  of 
this  thesis. 

Figures  B.1-B.3  show  the  initial  GRACE  residuals  for  each  component.  The 
error  growth  in  the  z  component  has  to  be  due  to  a  frequency  error  in  oq ,  since  u 2 
does  not  show  up  in  the  z  component.  The  error  growth  in  the  x  and  y  components 
may  be  due  to  the  error  in  oq,  or  an  error  in  uq,  or  possibly  both. 
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Figure  B.l:  The  residual  growth  over  time  in  the  x  component. 


Since  u;i  is  present  in  all  the  components,  this  frequency  was  changed  first. 
Examining  the  slope  of  the  error  growth  on  the  z  component,  the  author  found  the 
expected  error  to  be  on  the  order  of  10-5  rad/TU.  The  z  residuals  after  a  change  of 
IE  —  5  rad/TU  are  presented  in  Figure  B.4,  and  show  a  small  improvement  over  the 
initial  z  residuals  presented  in  Figure  B.3.  Adjustments  were  made  to  reduce  the  error 
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Figure  B.2:  The  residual  growth  over  time  in  the  y  component. 


15 


1 e+006 

800000 

600000 

400000 

200000 

0 

-200000 

-400000 

-600000 

-800000 

-1  e+006 


-15  -10  -5  0  5  10  15 


Days 

Figure  B.3:  The  residual  growth  over  time  in  the  z  component. 
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growth  in  z  over  time,  until  the  frequency  adjustment  reached  7E—5  rad/TU.  Further 
adjustments  on  the  order  of  10~5  rad/TU  resulted  in  the  error  growth  increasing 
again.  Figure  B.5  shows  the  residuals  after  the  frequency  adjustment. 
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Figure  B.4:  The  residual  growth  over  time  in  the  z  com¬ 
ponent  with  a  frequency  adjustment  of  8uj\  = 
lEnegb  rad/TU. 

Figure  B.5  shows  the  error  growth  in  x  and  y  decreased  from  those  shown  in 
Figures  B.l  and  B.2.  However,  the  residuals  still  show  a  frequency  error  on  the 
order  of  10-6  rad/TU.  Since  the  uj\  frequency  error  has  been  partially  resolved,  the 
remaining  frequency  error  in  x  and  y  must  be  from  errors  in  ui2-  Following  the  same 
method  as  above,  the  frequency  error  was  found  to  be  ~  5 E  —  6  rad/TU .  Figure  B.6 
shows  the  residuals  with  the  both  frequencies  adjusted. 

An  interesting  characteristic  shows  up  in  the  residuals  displayed  in  Figure  B.6. 
The  residuals  on  the  left  are  much  greater  than  the  residuals  on  the  right.  Finer 
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Figure  B.5:  The  GRACE  residuals  after  the  frequency  change 
of  5u\  =  7Eneg5  rad/TU. 


adjustments  to  the  frequencies  only  leads  to  balancing  the  ends  of  the  residuals, 
but  the  overall  magnitudes  are  not  greatly  reduced.  This  can  be  seen  in  Figure 
4.7.  The  reason  the  overall  residuals  cannot  be  furthered  improved  is  the  existence 
of  atmospheric  drag  on  the  GRACE  satellite.  This  also  accounts  for  the  quadratic 
shape  in  the  z  residuals. 
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Figure  B.6:  The  GRACE  residuals  after  the  frequency  changes 
of  Sui  =  7Eneg5  rad/TU  and  Su>2  = 
hEnegQ  rad/TU. 
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Appendix  C.  Additional  Jason- 1  Results 

This  appendix  provides  additional  results  for  the  Jason-1  satellite  that  were  obtained 
during  the  analysis.  A  brute  force  method  for  changing  the  x>2  frequency  is  highlighted 
here. 

Figures  C.1-C.3  show  the  difference  in  the  initial  x  and  y  residuals  and  the  initial 
2  residuals.  Notice  that  the  z  residuals  are  same  as  those  presented  in  Figure  4.15 
without  any  changes  to  the  frequencies,  but  the  x  and  y  are  off  significantly.  From 
this,  the  author  knew  that  the  only  significant  frequency  error  was  in  u/2- 
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Figure  C.l:  The  residual  growth  over  time  in  the  x  component. 


The  slope  of  the  residual  growth  in  x  and  y  indicates  that  the  error  in  u 2  is  on 
the  order  of  10"5.  Figure  C.4  shows  that  adjusting  the  frequency  by  IE  —  5  reduced 
the  residuals  by  about  However,  the  residuals  are  still  off  by  ~  100  km.  Figure 
C.5  shows  that  adjusting  the  frequency  by  2 E  —  5  reduced  the  residuals  by  an  order 
of  magnitude. 

The  goal  for  this  work  was  to  get  the  residuals  to  about  1  km.  Having  already 
improved  the  residuals  from  cs  200  km  down  to  ~  15  km,  the  author  assumed  that 
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Figure  C.2:  The  residual  growth  over  time  in  the  y  component. 
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Figure  C.3:  The  residual  growth  over  time  in  the  0  component. 
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Figure  C.4:  The  Jason-1  residuals  after  the  frequency  changes 
of  5u>2  =  lEnegh  rad/TU . 


further  frequency  improvements  would  be  of  lower  order.  The  first  attempt  was  to 
try  Slu2  =  2.1  E  —  5  rad/TU.  As  seen  in  Figure  C.6,  this  resulted  in  an  increase 
in  the  residuals.  This  meant  that  the  frequency  adjustment  was  actually  too  high. 
Continuing  to  adjust  the  frequencies  led  to  the  final  frequency  adjustment  discussed 
in  Section  4.2.3. 
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Figure  C.5:  The  Jason-1  residuals  after  the  frequency  changes 
of  Slu2  =  2Enegh  rad/TU. 
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Figure  C.6:  The  Jason-1  residuals  after  the  frequency  changes 
of  5u> 2  =  2.1Enegh  rad/TU. 
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